Regresja i analiza wariancji

Regresja liniowa i wieloraka

Matematyka stosowana

library(tidyverse)
library(ISLR)
library(kableExtra)
library(GGally)
library(knitr)
library(lmtest)
library(corrplot)
library(olsrr)
library(caret)
library(plotly)

Regresja liniowa

Przedstawienie oraz opis danych

house <- read.csv('C:/Users/Kuba/Downloads/house.csv', header = TRUE)
head(house)

Zbiór house przedstawia ceny domów w Stanach Zjednoczonych. Niektóre istotne kolumny w zbiorze danych:

• AvgAreaIncome - średni dochód powierzchniowy mieszkańców domu w [$],
• HouseAge - wiek domu,
• NumberOfRooms - liczba pokoi (średnia),
• NumberOfBedooms - liczba sypialni (średnia),
• AreaPopulation - liczba ludności w dzielncy/regionie, w której znajduje się dany dom,
• Price - cena domu w [$].

kable(summary(house), escape = F) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

	Avg..Area.Income	House.Age	Number.of.Rooms	Number.of.Bedrooms	Area.Population	Price	Address
	Min. : 17797	Min. :2.644	Min. : 3.236	Min. :2.000	Min. : 172.6	Min. : 15939	Length:4548
	1st Qu.: 61485	1st Qu.:5.332	1st Qu.: 6.300	1st Qu.:3.140	1st Qu.:29423.2	1st Qu.: 997775	Class :character
	Median : 68817	Median :5.961	Median : 7.002	Median :4.050	Median :36215.6	Median :1234571	Mode :character
	Mean : 68612	Mean :5.979	Mean : 6.988	Mean :3.982	Mean :36187.5	Mean :1233916	NA
	3rd Qu.: 75821	3rd Qu.:6.658	3rd Qu.: 7.666	3rd Qu.:4.490	3rd Qu.:42880.6	3rd Qu.:1470616	NA
	Max. :107702	Max. :9.519	Max. :10.760	Max. :6.500	Max. :69592.0	Max. :2469066	NA

Założenia

Zależność liniowa

Przekształcamy “tibbl’a” ze wszystkimi zmiennymi na mniejszą wersję, zawierającą wyłącznie zmienne ciągłe (usuwamy kolumnę Address).

data(house)

## Warning in data(house): data set 'house' not found

ggpairs(house[, c(1:6)])

Powyższy wykres przedstawia zależności pomiędzy wszystkimi ciągłymi zmiennymi. Można zauważyć, że istnieje kilka znacznych korelacji. Jedna, która wyróżnia się ponad innymi to ta pomiędzy zmienną Price i AvgAreaIncome ze współczynnikiem korelacji Pearsona równym \(0,641\). Korelacje te są dość istotne, ponieważ wartość współczynnika jest bliska \(1\) dla współczynnika dodatniego. Dlatego ceny domów w Stanach Zjednoczonych będziemy przewidywać na podstawie zmiennej Avg..Area.Income. Innymi słowy zmienna Price to zmienna objaśniana, a Avg..Area.Income objaśniająca.

Rozkłady zmiennych

Sprawdzimy teraz rozkład tych dwóch zmiennych za pomocą histogramów.

ggplot(house, aes(x=Price))+
  labs(title = 'Histogram cen domów',
       x = 'Cena',
       y = 'Częstotliwość') + 
  geom_histogram(aes(y=..density..),
                 color="khaki1", 
                 fill="mediumslateblue", 
                 binwidth = 50000) +
  geom_density(alpha=.45, fill="springgreen4") +
  theme_classic()+
  theme(plot.title = element_text(hjust = 0.5))

## Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
## i Please use `after_stat(density)` instead.

ggplot(house, aes(x=Avg..Area.Income))+
  labs(title = 'Histogram średniego dochodu powierzchniowego w domu',
       x = 'Średni dochód powierzchniowy w domu',
       y = 'Częstotliwość') + 
  geom_histogram(aes(y=..density..),
                 color="khaki1", 
                 fill="mediumslateblue", 
                 binwidth = 2100) +
  geom_density(alpha=.45, fill="springgreen4") +
  theme_classic()+
  theme(plot.title = element_text(hjust = 0.5))

Wizualizując rozkłady tych dwóch zmiennych możemy domyśleć się, że posiadają rozkład normalny. Dla potwierdzenia naszej hipotezy użyjemy testu Shapiro-Wilka, który zakłada, że zmienna posiada rozkład normalny.

shapiro.test(house$Price)

## 
##  Shapiro-Wilk normality test
## 
## data:  house$Price
## W = 0.9998, p-value = 0.9688

shapiro.test(house$Avg..Area.Income)

## 
##  Shapiro-Wilk normality test
## 
## data:  house$Avg..Area.Income
## W = 0.99958, p-value = 0.4456

Wnioskując, w obu przypadkach otrzymaliśmy \(p\)-value \(> 0,05\), zatem nie mamy podstaw do odrzucenia \(H_0\).

scatter_plot <- plot_ly(data = house, x = ~Avg..Area.Income, y = ~Price, color = ~Price) %>%
layout(title = 'Wykres punktowy zależności ceny domu od średniego dochodu powierzchniowego', plot_bgcolor = "#e5ecf6")
scatter_plot 

## No trace type specified:
##   Based on info supplied, a 'scatter' trace seems appropriate.
##   Read more about this trace type -> https://plotly.com/r/reference/#scatter

## No scatter mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode

Wykres nie wyklucza zależności liniowej pomiędzy zmiennymi Price i Avg..Area.Income, więc przejdziemy do stworzenia modelu regresji liniowej tych zmiennych.

Model regresji liniowej

lmfit1 <- lm(Price~Avg..Area.Income, data = house)
summary(lmfit1)

## 
## Call:
## lm(formula = Price ~ Avg..Area.Income, data = house)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -965631 -187471   -3066  184157  982367 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -2.262e+05  2.622e+04  -8.627   <2e-16 ***
## Avg..Area.Income  2.128e+01  3.775e-01  56.365   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 272100 on 4546 degrees of freedom
## Multiple R-squared:  0.4114, Adjusted R-squared:  0.4112 
## F-statistic:  3177 on 1 and 4546 DF,  p-value: < 2.2e-16

Z naszego modelu, możemy wywnioskować:

• współczynnik kierunkowy \(\beta_1\) w naszym modelu regresji wynosi \(21,28\), natomiast \(\beta_0\) czyli wyraz wolny \(-226200\), czyli wraz ze wzrostem średniego dochodu powierzchniowego w domu rośnie jego cena,
• \([***]\) przy współczynniku \(\beta_1\) oznacza, że średni dochód powierzchniowy w domu ma istotny wpływ na cenę,
• błąd resztowy (odchylenie standardowe składnika resztowego) wynosi \(272100 \; \$\) co oznacza, że wartości obliczone na podstawie modelu różnią się od rzeczywistości średnio \(\pm\; 272100 \; \$\),
• współczynnik determinacji \((R^2)\) informuje nas o tym, jaki procent wariancji zmiennej objaśnianej został wyjaśniony przez funkcję regresji, Współczynnik \(R^2\) (multiple R-squared) wynosi \(0,4114\) czyli średni dochód powierzchniowy mieszkańców domu wyjaśnia \(41,14\%\) zmienności ceny.

Rozkład reszt

ggplot(lmfit1, aes(x=resid(lmfit1)))+
  labs(title = 'Histogram reszt modelu',
       x = 'Reszty',
       y = 'Częstotliwość') + 
  geom_histogram(aes(y=..density..),
                 color="khaki1", 
                 fill="mediumslateblue", 
                 binwidth = 50000) +
  geom_density(alpha=.45, fill="springgreen4") +
  theme_classic()+
  theme(plot.title = element_text(hjust = 0.5))

ggplot(lmfit1, aes(sample=lmfit1$residuals)) + 
  geom_qq() + 
  geom_qq_line(color = 'red') + 
  labs(title='Wykres kwartyl-kwartyl reszt', x='Kwartyle teoretyczne', y='Kwartyle próbkowe')+
  theme(plot.title = element_text(hjust = 0.5))

shapiro.test(lmfit1$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  lmfit1$residuals
## W = 0.99971, p-value = 0.7937

Wszystkie metody użyte przez nas świadczą o tym, że reszty naszego modelu posiadają rozkład normalny.

Zerowa średnia reszt

Do sprawdzenia zerowej średniej reszt użyjemy testu t-studenta.

t.test(lmfit1$residuals)

## 
##  One Sample t-test
## 
## data:  lmfit1$residuals
## t = -9.1167e-15, df = 4547, p-value = 1
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -7908.126  7908.126
## sample estimates:
##     mean of x 
## -3.677457e-11

W tym wypadku test t-student wykazał, że średnia jest równa zero.

ggplot(lmfit1, aes(.fitted, .resid)) + 
  geom_point() + 
  stat_smooth(method='loess', formula=y~x, se=F) +
  geom_hline(yintercept=0, linetype='dashed', color='red') +
  labs(title='Wykres zależności reszt od dopasowanych wartości', 
       x='Dopasowane wartości',
       y='Reszty')+
  theme(plot.title = element_text(hjust = 0.5))

Patrząc na wykres zależności reszt od dopasowanych wartości możemy zauważyć większe odchylenia dla małych jak i dużych wartości zmiennej objaśnianej (Price).

Niezależność reszt

W celu sprawdzenia niezależności reszt przeprowadzimy test Durbin-Watson’a.

lmtest::dwtest(lmfit1)

## 
##  Durbin-Watson test
## 
## data:  lmfit1
## DW = 1.9676, p-value = 0.1373
## alternative hypothesis: true autocorrelation is greater than 0

W naszym przypadku \(p\)-value \(> 0,05\), więc nie mamy podstaw, aby odrzucić hipotezę o niezależności w resztach.

Homoskedastyczność

Aby sprawdzić homoskedastyczność posłuży nam wykres zależności pierwiastka standaryzowanych reszt od dopasowanych wartości.

lmfit1 %>%
ggplot(aes(.fitted, sqrt(abs(.stdresid)))) + 
  geom_point() + stat_smooth(method='loess', formula=y~x, se=F) +
  labs(title='Zależność pierwiastka standaryzowanych reszt od dopasowanych wartości', x='Dopasowane wartości', y='Pierwiastek standaryzowanych reszt')+
  theme(plot.title = element_text(hjust = 0.5))

Rozrzut reszt na wykresie jest mniej więcej równy dla wszystkich dopasowanych wartości. Widzimy odchylenie na początku niebieskiej prostej. Jest ono spowodowane wartością odstającą. Deszty są równomiernie rozrzucone wokół niebieskiej linii. Dla pewności przeprowadzimy test Breusch-Pagan’a, który za hipotezę zerową zakłada homoskedastyczność reszt.

lmtest::bptest(lmfit1)

## 
##  studentized Breusch-Pagan test
## 
## data:  lmfit1
## BP = 1.893, df = 1, p-value = 0.1689

\(P\)-value \(> 0,05\), więc zakładamy homoskedastyczność reszt.

Podsumowanie

Słowem podsumowania stwierdzamy, że powyższy model regresji liniowej jest zgodny z naszymi założeniami. Korzystająć ze współczynnika korelacji Pearsona dopasowaliśmy zmienną objaśniającą (Avg..Area.Income). Przewidywana zmienna Price jest najbardziej zależna od zmiennej Avg..Area.Income. Analiza reszt pokazała, że model jest zgodny założeniami i ma postać:

\(Price = -226200 + 21,28 \cdot Avg..Area.Income\)

---

Regresja wieloraka

Spójrzmy raz jeszcze na nasz zbiór danych.

head(house)

Ze względnu na to, że w naszym modelu nie może byc zmiennych jakościowych, wybieramy tylko zmienne ilościowe przy pomocy funkcji select().

new_house <- house %>% 
                select(Avg..Area.Income, House.Age, Number.of.Rooms, Number.of.Bedrooms, Area.Population, Price)

head(new_house)

Model regresji wielorakiej

lmfit2 <- lm(Price~., data = new_house)
summary(lmfit2)

## 
## Call:
## lm(formula = Price ~ ., data = new_house)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -338141  -70192    -101   69104  347218 
## 
## Coefficients:
##                      Estimate Std. Error  t value Pr(>|t|)    
## (Intercept)        -2.639e+06  1.791e+04 -147.400   <2e-16 ***
## Avg..Area.Income    2.161e+01  1.404e-01  153.929   <2e-16 ***
## House.Age           1.653e+05  1.514e+03  109.228   <2e-16 ***
## Number.of.Rooms     1.211e+05  1.682e+03   72.003   <2e-16 ***
## Number.of.Bedrooms  1.458e+03  1.377e+03    1.059    0.289    
## Area.Population     1.519e+01  1.513e-01  100.357   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 101100 on 4542 degrees of freedom
## Multiple R-squared:  0.9188, Adjusted R-squared:  0.9187 
## F-statistic: 1.028e+04 on 5 and 4542 DF,  p-value: < 2.2e-16

Dzięki funkcji summary() odczytujemy, że \(R^2\) naszego modelu wynosi \(0,9187\). Tak zbudowany model wyjaśnia \(91,87\%\) zmienności bieżącej ceny, ale nie wszystkie zmienne są w tym modelu istotne.

Zinterpretujmy nasze wyniki:

• średni dochód mieszkańców domu zwiększa jego cenę średnio o \(21,61 \; \$\),
• im mniejszy wiek domu tym zwiększa się jego wartość średnio o \(165300 \; \$\),
• ilość pokoi w domu zwiększa jego wartość średnio o \(121100 \; \$\),
• ilość sypialni w domu zwiększa jego wartość średnio o \(1458 \; \$\),
• liczność populacji zwiększa wartość domu średnio o \(15,19 \; \$\),

Algorytm AIC

Wykorzystamy współczynnik informacyjny Akaike'a (AIC). Jest to estymator błędu przewidywania. Biorąc pod uwagę zbiór modeli dla danych, AIC szacuje jakość każdego modelu, w stosunku do każdego z pozostałych modeli. W ten sposób AIC dostarcza środków do wyboru modelu. Dzięki niemu uzyskamy najlepszą możliwą dla nas jakość modelu poprzez eliminację następnych zmiennych. Skorzystamy z powyższego algorytmu przy pomocy funkcji step(). Algorytm będzie działał na zasadzie procedury eliminacji wstecznej.

 step(lmfit2, direction = "backward")

## Start:  AIC=104827.4
## Price ~ Avg..Area.Income + House.Age + Number.of.Rooms + Number.of.Bedrooms + 
##     Area.Population
## 
##                      Df  Sum of Sq        RSS    AIC
## - Number.of.Bedrooms  1 1.1473e+10 4.6440e+13 104827
## <none>                             4.6429e+13 104827
## - Number.of.Rooms     1 5.2996e+13 9.9425e+13 108289
## - Area.Population     1 1.0295e+14 1.4938e+14 110140
## - House.Age           1 1.2196e+14 1.6839e+14 110685
## - Avg..Area.Income    1 2.4220e+14 2.8863e+14 113136
## 
## Step:  AIC=104826.6
## Price ~ Avg..Area.Income + House.Age + Number.of.Rooms + Area.Population
## 
##                    Df  Sum of Sq        RSS    AIC
## <none>                           4.6440e+13 104827
## - Number.of.Rooms   1 6.8538e+13 1.1498e+14 108948
## - Area.Population   1 1.0295e+14 1.4939e+14 110138
## - House.Age         1 1.2204e+14 1.6848e+14 110685
## - Avg..Area.Income  1 2.4255e+14 2.8899e+14 113139

## 
## Call:
## lm(formula = Price ~ Avg..Area.Income + House.Age + Number.of.Rooms + 
##     Area.Population, data = new_house)
## 
## Coefficients:
##      (Intercept)  Avg..Area.Income         House.Age   Number.of.Rooms  
##       -2.640e+06         2.161e+01         1.654e+05         1.220e+05  
##  Area.Population  
##        1.519e+01

Jak widzimy kolejno usuwamy zmienne, które w naszym modelu mają najmniejszą wartość AIC, najpierw był to Number.of.Bedrooms. Dochodzimy do momentu w algorytmie gdzie najniższą wartość AIC ma <none> zatem poniższe zmienne zostają w naszym modelu regresji wielorakiej. Widzimy, że dla wszystkich zmiennych \(R^2\) wynosi \(0,9187\). Sprawdźmy teraz ile wynosi nasz \(R^2\) dla naszego gotowego modelu, po redukcji tej zmiennej zmiennych:

new_lmfit2 <- lm(Price~ Avg..Area.Income + House.Age + Area.Population + Number.of.Rooms, data= new_house)
summary(new_lmfit2)

## 
## Call:
## lm(formula = Price ~ Avg..Area.Income + House.Age + Area.Population + 
##     Number.of.Rooms, data = new_house)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -339350  -70321    -126   69246  347378 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -2.640e+06  1.790e+04 -147.45   <2e-16 ***
## Avg..Area.Income  2.161e+01  1.403e-01  154.04   <2e-16 ***
## House.Age         1.654e+05  1.513e+03  109.27   <2e-16 ***
## Area.Population   1.519e+01  1.513e-01  100.35   <2e-16 ***
## Number.of.Rooms   1.220e+05  1.490e+03   81.88   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 101100 on 4543 degrees of freedom
## Multiple R-squared:  0.9188, Adjusted R-squared:  0.9187 
## F-statistic: 1.284e+04 on 4 and 4543 DF,  p-value: < 2.2e-16

Nasz współczynnik \(R^2\) dla obu modeli jest cały czas taki sam, zatem zmienna Number.of.Bedrooms nie ma wpływu istotnie na jakość modelu. Natomiast pozostałe zmienne są istotne dlatego nie odrzuciliśmy ich przy użyciu algorytmuAIC.

Przykład z usunięciem zmiennej objaśniającej

new_lmfit2_no_avg <- lm(Price~ House.Age  + Area.Population  + Number.of.Rooms, data=new_house)
summary(new_lmfit2_no_avg)

## 
## Call:
## lm(formula = Price ~ House.Age + Area.Population + Number.of.Rooms, 
##     data = new_house)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -895558 -167437    2368  171008  878807 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -1.126e+06  3.732e+04  -30.16   <2e-16 ***
## House.Age        1.630e+05  3.775e+03   43.18   <2e-16 ***
## Area.Population  1.501e+01  3.774e-01   39.77   <2e-16 ***
## Number.of.Rooms  1.205e+05  3.715e+03   32.42   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 252200 on 4544 degrees of freedom
## Multiple R-squared:  0.4945, Adjusted R-squared:  0.4941 
## F-statistic:  1481 on 3 and 4544 DF,  p-value: < 2.2e-16

Jak widzimy powyżej, nasz \(R^2\) znacznie zmalał, co oznacza iż ten nowy model objaśnia tylko \(49,41\%\) zmienności bieżacej ceny domów. Jest to niekorzystne w wypadku gdy możemy objaśniać aż \(91,87\%\) ceny.

Założenia

Badanie współliniowości

Omówmy temat współliniowości zmiennych objaśniających:

ols_vif_tol(new_lmfit2)

Kolumna Tolerance wskazuje wartość procentową niewyjaśnionej zmienności danej zmiennej przez pozostałe zmienne objaśniające. Dla przykładu: współczynnik tolerancji dla Avg..Area.Income wynosi $0.9997982 $ co oznacza, że \(99,97\%\) zmienności Avg..Area.Income nie jest wyjaśnione przez pozostałe zmienne w modelu. Kolumna VIF jest obliczana na podstawie wartości współczynnika tolerancji i pokazuje o ile wariancja szacowanego współczynnika regresji jest podwyższona z powodu współliniowości danej zmiennej objaśniającej z pozostałymi zmiennymi objaśniającymi. Wszystkie zmienne mają VIF na podobnym poziomie w granicach \(1.0\) zatem nie wskazują współliniowości (VIF \(> 4\) wskazuje na współliniowość zmienncyh, także u nas takie zmienne nie występują).

Normalność reszt

Sprawdźmy normalność naszych reszt:

ggplot(lmfit2, aes(x=resid(lmfit2)))+
  labs(title = 'Histogram reszt modelu',
       x = 'Reszty',
       y = 'Częstotliwość') + 
  geom_histogram(aes(y=..density..),
                 color="khaki1", 
                 fill="mediumslateblue", 
                 binwidth = 18000) +
  geom_density(alpha=.2, fill="springgreen4") +
  theme_classic()+
  theme(plot.title = element_text(hjust = 0.5))

Reszty w naszym modelu wydają się być zbliżone do rozkładu normalnego, natomiast nie możemy tego stwierdzić w 100%. Sprawdźmy wykres QQ-plot do zbadania rozkładu normalnego naszych zmiennych.

ggplot(new_lmfit2, aes(sample=new_lmfit2$residuals)) + 
  geom_qq() + 
  geom_qq_line(color = 'red') + 
  labs(title='Wykres kwartyl-kwartyl reszt', x='Kwartyle teoretyczne', y='Kwartyle próbkowe')+
  theme(plot.title = element_text(hjust = 0.5))

Widzimy, że prawie wszystkie punkty leżą na czerwonej lini, jedynie problem jest na początku oraz na końcu.

Zależność reszt

new_lmfit2 %>%
ggplot(aes(.fitted, .resid)) + 
  geom_point() + 
  stat_smooth(method='loess', formula=y~x, se=F) +
  geom_hline(yintercept=0, linetype='dashed', color='red') +
  labs(title='Wykres zależności reszt od dopasowanych wartości', 
       x='Dopasowane wartości',
       y='Reszty')+
  theme(plot.title = element_text(hjust = 0.5))

Na wykresie zależności reszt od dopasowanych wartości możemy zauważyć delikatne odchylenie w początkowej jego fazie. Innymi słowy, współczynniki modelu regresji powinny być godne zaufania i nie musimy wykonywać transformacji na danych.

Identyfikacja wartości odstających

ols_plot_cooksd_bar(new_lmfit2)

ols_plot_resid_stud_fit(new_lmfit2)

Przypisanie tej funkcji do obiektu zwraca nam tabelę z numerami zidentyfikowanych obserwacji wpływowych. Sposób ten pozwala nam na identyfikację punktów, które negatywnie wpływają na model regresji są one oznaczone kolorem czerwonym. Miara jest kombinacją wartości dźwigni i reszt każdej obserwacji; im wyższa dźwignia i reszty, tym wyższa odległość Cooka.

Homoskedastyczność

new_lmfit2 %>%
ggplot(aes(.fitted, sqrt(abs(.stdresid)))) + 
  geom_point() + stat_smooth(method='loess', formula=y~x, se=F) +
  labs(title='Zależność pierwiastka standaryzowanych reszt od dopasowanych wartości', x='Dopasowane wartości', y='Pierwiastek standaryzowanych reszt')

Na powyższym wykresie zależności pierwiastka standaryzowanych reszt od dopasowanych wartości, możemy zauważyć, iż wartości są rozmieszczone mniej więcej w takim samym odstępie od niebieskiej linii. Nie wykluczamy, więc homoskedastyczności. Dla pewności przeprowadźmy test Breusch-Pagan’a:

lmtest::bptest(new_lmfit2)

## 
##  studentized Breusch-Pagan test
## 
## data:  new_lmfit2
## BP = 5.8309, df = 4, p-value = 0.2121

Nie mamy wystarczających dowodów, aby stwierdzić, że heteroskedastyczność występuje w naszym modelu regresji, zatem możemy założyć homoskedastyczność.

Podsumowanie

Stworzony przez nas model spełnia wszystkie założenia. Korzystając z algorytmu AIC dobraliśmy najlepszy możliwy model regresji wielorakiej, który dany jest wzorem:

\(Price = -2640000 + 21,61 \cdot Avg..Area.Income + 165400 \cdot House.Age +\\+ 15,19 \cdot Area.Population + 122000 \cdot Number.of.Rooms\)